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1  Background 


The  reliable,  accurate,  and  efficient  computational  modeling  of  plasma  dynamics  remains  a  challenge 
of  very  significant  proportions  but  also  with  a  broad  range  of  possible  applications,  e.g.,  high-power 
microwave  generation;  large  scale  particle  accelerators;  fusion  energy,  both  by  means  of  magnetic 
confinement  and  laser  ignited  devices;  and  a  variety  of  plasma  based  technology  as  well  as  numerous 
problems  of  mainly  military  interest.  This  warrants  that  significant  resources  be  spend  on  the 
continued  development  of  accurate,  robust,  and  efficient  tools  for  the  modeling  of  this  broad  range 
of  phenomena  and  scales. 

The  direct  and  straightforward  solution  of  many  of  the  problems,  at  least  in  the  collisionless 
region,  could  in  principle  be  accomplished  by  solving  the  phase-space  Vlasov  equation 

^+v-Vxf  +  F-Vvf  =  0  ,  (1) 

where  f(x,v,t)  is  the  probability  density,  one  for  each  charged  particle  species,  and  the  force,  F, 
is  the  electromagnetic  Lorenz  force 

F  =  —[E  +  vxB]  , 
m 

for  each  particle  type  of  charge  q  and  mass  m.  In  case  of  relativistic  speeds,  the  mass  needs  to 
be  corrected  according  to  special  relativity.  The  electromagnetic  fields,  ( E,B )  satisfy  Maxwell 
equations 


dfdxEt  —  Vx  x  B  +  J  ,  — -  =  ~VX  x  E  ,  (2) 

at 

where  we  have  assumed  vacuum  conditions  at  unity  for  simplicity,  i.e.,  e0  =  l,/j,0  =  1,  and  the  speed 
of  light  is  normalized  to  one.  The  fields  are  further  constrained  through 


V*  •  E  =  p  ,  Vj.  •  B  =  0  , 

where  the  former  ensures  charge  conservation  through 


(3) 


dp 


•  J  =  0 


for  the  space  charge  p(x,t)  and  current  density  J(x,f). 

This  system  needs  to  be  solved  in  the  fully  dynamic  range  for  cases  where  the  particles  are  close 
to  relativistic  speeds  as  is  the  case  for  many  areas  of  modern  technology,  e.g.,  high-power  microwave 
generation,  accelerator  modeling  and  fusion  applications. 

However,  the  6+1  dimensional  nature  of  Eq.(l)  and  the  very  complex  phase  space  dynamics, 
essentially  implies  that  direct  modeling  remains  out  of  reach  for  problems  in  complex  geometries 
and  of  realistic  complexity.  On  the  other  hand,  this  approach  offers  a  direct  measure  stick  against 
which  to  compare  new  methods  and  evaluate  their  performance. 

The  practical  difficulties  associated  with  the  direct  solution  of  the  above  system  has  encouraged 
the  development  of  a  number  of  indirect  approaches  to  solve  this  type  of  problems.  The  most  widely 
used  approach  is  the  particle-in-cell  method  (PIC)  which  is  based  on  a  sampling  technique.  Here, 
one  samples,  using  a  large  number  of  particles,  the  local  distribution  function,  f(x,v,t),  in  velocity 
space  and  then  move  these  particles  around  in  a  fully  Lagrangian  fashion  using  Newton’s  law  of 
motion.  In  other  words,  one  represents  the  first  two  moments,  p{x,t),  and  J(x,t),  as 
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p(x,t)  =  J2Pkix,t)  ’  J(abf)  =  H*j*(abf)  > 

k  k 

where  pk  and  Jj.  represent  the  charge  and  and  current  originating  from  particle  k.  Each  of  these  k 
particles  are  them  moved  according  to 

=pk[E(xk)  +vk  x  B(xfc)]  . 

Thus,  the  computational  work  is  shifted  from  a  6+1  dimensional  problem  to  a  3+1  problem  but 
with  a  large  number  of  particles  moving  freely  around  and  for  which  there  is  constant  interaction 
between  the  fields  and  the  particles.  The  simple  sampling  approach  taken  in  this  technique  also 
means  that  a  very  large  number  of  particles  has  to  be  used  to  limit  the  classic  1/K  sampling  noise. 

While  these  methods  generally  are  efficient  and  has  served  the  community  well,  they  typically 
rely  on  the  use  of  relatively  simple  grids  and  low  order  accurate  methods,  in  particular  for  the  fully 
electromagnetic  problem,  for  the  field  solvers  -  typically  2nd  order  finite  difference  or  finite  element 
methods,  the  latter  mainly  for  the  low-speed  case.  As  a  consequence  of  this,  many  of  these  algorithms 
are  now  being  stretched  due  to  the  quest  to  resolve  problems  with  a  large  range  of  active  scales, 
both  in  time  and  space.  Further  complications  arise  when  the  interaction  between  geometries  and 
the  fields  become  important  such  as  in  magnetrons  and  klystons  for  microwave  generation,  as  well 
as  problems  associated  with  inherent  numerical  properties  of  the  low  order  schemes,  e.g,  numerical 
dispersion  and  numerical  Cherenkov  radiation. 

In  this  proposal  we  have  consider  the  continued  development  and  maturation  of  a  computational 
kernel  based  on  the  discontinuous  Galerkin  method  (DG-FEM)  for  the  field  solver  while  various  al¬ 
ternatives  are  possible  for  the  particle  models.  This  approach  has  several  advantages  for  applications 
of  the  type  discussed  in  this  effort,  i.e. ,  their  ability  to  easily  handle  complicated,  electrically  large 
geometries  and  complex  boundary  conditions,  their  flexibility  to  support  nonconforming  hp  adap¬ 
tivity,  and  the  ability  to  efficiently  deal  with  both  the  high-speed  fully  electromagnetic  problem  and 
the  low  speed  problem  within  one  formulation. 

DG-FEM  have  experienced  an  explosion  in  interest  over  the  last  decade,  although  with  relatively 
limited  impact  in  plasma  physics  applications.  This  is  in  particular  the  situation  for  particle-mesh 
methods  where,  as  discussed  above,  classic  2nd  order  accurate  PIC  methods  remain  the  dominating 
tool  in  spite  of  well  known  inherent  problems. 

As  we  shall  discuss  shortly,  we  have  made  significant  progress  towards  the  development  of  an 
alternative  to  the  current  simulation  tools.  Ultimately,  the  hope  is  that  the  continuation  of  this 
effort  will  result  in  a  software  tool  of  sufficient  robustness  and  versatility  to  be  useful  to  AFRL 
researchers  in  the  development  and  evaluation  of  new  technology,  critical  to  the  mission  of  USAF. 


2  Main  Developments 

A  thoroughly  validated,  robust,  and  efficiently  solver  for  the  electromagnetics  problem  continues 
to  be  extended  to  provide  a  basic  field  solver  for  the  EM-PIC  problem.  This  has  resulted  in  a 
first  generation  of  a  DG-FEM  based  PIC  method  on  unstructured  grids.  The  particle  mover  relies 
on  high-order  interpolation,  efficient  local  search  algorithms  to  locate  the  particles,  and  a  level  set 
approach  to  implement  elastic/inelastic  interactions  with  boundaries.  Charge  and  current  redis¬ 
tribution  computations  use  large  smooth  weight  functions,  dramatically  decreasing  the  number  of 
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particles  needed  in  a  computational  cell  as  compared  to  the  simple  redistribution  schemes  typically 
used.  A  large  smooth  Gaussian-like  distribution  function  is  used  as  a  regularized  for  the  particle 
per  particle  to  sums  the  influence  of  the  particles  into  the  charge  and  current  density  field  variables. 
Currently,  a  constant  particle  influence-radius  ensures  charge  conservation  in  the  domain,  but  re¬ 
quires  a  variable  number  of  elements  per  particle  for  weighing  purposes  depending  on  the  size  of  the 
elements.  The  element's  that  are  being  influenced  by  a  particle  are  determined  in  a  pre-processing 
stage  by  finding  the  elements  connected  to  the  node  nearest  to  the  particle. 

In  this  effort,  we  have  primarily  concerned  ourselves  with  three  topics 

•  Time- advancement.  We  have,  with  great  success,  implemented  a  implicit-explicit  Runge-Kutta 
methods  [?]  in  which  the  full  EM  part,  Eq.(2),  is  done  implicitly,  while  the  charge  dynamics 
and  current  deposition  is  done  explicitly.  This  removes  all  grid  induced  stiffness  and  enables 
one  to  effectively  model  geometrically  complex  problems.  However,  stiffness  introduced  by 
physical  properties,  i.e.,  high  electron  cyclotron  frequencies  or  small  Debye  shielding  lengths 
still  imposes  certain  constraints  on  the  time  step.  This  has  lead  to  significant  speedups  for 
classes  of  problems  where  hyperbolic  divergence  cleaning  is  the  preferred  approach.  Further¬ 
more,  this  approach  is  well  suited  for  problems  in  complex  geometries. 

To  completely  overcome  stiffness,  including  that  introduced  by  physical  phenomena,  we  have 
also  considered  a  full  implicit  moment  method.  This  is  highly  efficient  for  problems  in  strong 
magnetic  fields  and  for  problems  with  particles  of  different  charge  and  mass,  e.g.,  electron-ion 
mixtures.  Currently,  this  method  is  at  best  2n  dorder  accurate  and  we  are  working  to  address 
this 

Both  methods  have  been  carefully  validated  and  their  design  accuracy  confirmed  for  simple 
test  problems,  as  exemplified  in  Fig  1. 

•  Techniques  for  phase-spase  control  A  key  concern  in  methods  of  this  type  is  the  very  high 
number  of  particles  typically  used.  This  normally  controls  the  computational  cost  entirely. 
We  have  pursue  a  solution  to  this  by  investigating  alternative  ways  of  representing  particles, 
techniques  for  splitting  and  combining  particles  and,  finally,  ways  of  estimating  the  phase- 
space  error. 

The  most  fruitful  approach  appears  to  be  that  of  solving  a  problem  twice  in  different  ways 
and  then  comparing  they  two.  For  this,  one  needs  an  alternative  equation  and  we  propose  to 
consider  the  <5/-equation 

85  f 

-gj: -  +  «  ■  Vx5f  +  F  ■  Wv6f  =  -V  ■  Vx/o  +  F  •  V„/0  , 

i.e.,  it  is  the  Vlasov  equation  (locally)  linearized  around  some  fo-  Clearly,  if  fo  is  kept  fixed 
in  time,  this  will  be  expensive  to  solve.  However,  if  we  use  non-parametric  estimation  to  re¬ 
evaluate  fo  at  fixed  temporal  intervals,  one  can  solve  the  above  linearized  Vlasov  equation  with 
a  low  number  of  particles.  The  predicted  changes  in  5f  will  provide  a  measure  of  the  phase- 
space  dynamics/sensitivity  and  we  can  then  compare  this  with  the  current  local  sampling  of 
phase-space  to  recover  a  measure  of  whether  phase-space  is  adequately  sampled. 

An  attractive  alternative  to  this  is  to  explore  appropriate  fluid  models  for  the  systems  being 
modeled  and  exploit  the  fact  that  DG-FEM  is  ideally  suited  for  solving  such  fluid  models, 
often  given  as  conservation  laws.  Furthermore,  since  the  cost  of  solving  the  fluid  model  is 
much  less  that  the  kinetic  model,  we  propose  to  run  both  models  simultaneously.  The  output 
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—  Runge-Kutta 

—  Implicit  moment 


Figure  1:  Direct  comparison  of  a  DG-FEM  based  scheme  with  two  different  implicit  temporal 
integration  schemes. 


from  the  fluid  model  then  provides  information  about  the  first  few  moments  of  the  kinetic 
simulation  and  a  comparison  between  the  results  from  the  fluid  model  and  the  kinetic  model 
can  be  used  to  identity  regions  close  to  equilibrium,  hence  requiring  just  a  few  particles,  and 
strongly  kinetic  regions  where  more  particles  are  needed.  Since  the  fluid  solver  is  used  for  error 
estimation  only,  it  can  be  solved  with  marginal  resolution  and  restarted  with  initial  conditions 
from  the  kinetic  solver  at  given  interval.  This  will  also  provide  a  suitable  return  place  in  case 
the  estimator  determines  that  the  physics  is  severely  under  resolved. 

Both  of  these  techniques  are  currently  being  explored  further  in  a  carefully  designed  experi¬ 
ment  to  judge  which  of  the  two  is  advantageous. 

•  Thorough  validation  The  scheme  has  been  implemented  and  extensively  tested  in  two  spatial 
dimensions,  both  for  low-speed  and  high-speed  test  cases  and  applications,  including  classic 
cases  such  as  plasma  waves,  Landau  damping,  two-stream  instabilities,  relativistic  Weibel 
instabilities  etc.  In  all  these  cases,  careful  comparisons  have  been  done  with  standard  methods 
based  on  finite  difference  methods,  yields  very  good  agreement.  An  example  is  shown  in  Fig. 
2  where  we  show  the  direct  comparison  for  the  relativistic  Weibel  instability 

Following  this  extensive  testing,  we  have  demonstrated  the  flexibility  and  robustness  of  the 
scheme  for  two-dimensional  simplification  closer  to  problems  of  realistic  complexity.  An  ex¬ 
ample  of  this  is  shown  in  Fig.3  where  we  show  the  modeling  of  a  relativistic  magnetron. 

Many  other  benchmarks  have  been  considered,  including  bore-magnetrons  and  non- dissipative 
magnetic  recognection. 
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-  X=2,  N=64,  Np=500 

.  X=10,  N=64,  Np=500 

-  X=2,  N=128,  Np=1000 

.  X=10,N=128,Np=1000 

- FDTD,  N=256 


Figure  2:  Direct  comparison  of  a  DG-FEM  based  scheme  with  a  classic  finite  difference  based  scheme 
for  modeling  of  the  Weibel  instability. 
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Figure  3:  Left)  Grid  for  2D  version  of  A6  relativistic  magnetron.  Right)  Charge  density  in  the 
magnetron,  running  in  the  2tt  mode.  Both  results  computed  using  a  DG-FEM  for  full  PIC 


G.  Jacobs,  G.  Lapenta,  and  J.S.  Hesthaven,  2006,  High-Order  Nodal  Discontinuous  Galerkin 
Particle-in- Cell  Methods  for  Modeling  of  Weibel  Instabilities,  AIAA  paper,  Reno,  2006. 

Two  more  papers  are  currently  being  finished,  devoted  to  the  development  and  validation 
of  the  time-stepping  schemes,  and  a  detailed  validation  study  of  geomagnetic  reconnection 
benchmark  (GEM). 
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